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Field of the Invention 

ring of digitally represented samples of analog data can 
,erfonmng the required processing in some transformed 
1 and its inverse to simplify the distal or other manip- 
orms the data from its initial form to an intermediate 
anipulation, and finally appKes the inverse transform to 
rm from its intermediate form. By way of example and 
nitial form of the data may be digital sanjples of anar 
form of the data may be the coefficients of the Fourier 
Le data may be the coefficients of the Fourier transform, 
le data may be the coefficients in some wavelet or other 
he data may be an analog signal, and the interme.hate 
diotal representation of it; the initial form of the data 
lalog signal, and the intermediate form of the fata may 
1 samples. For instance, transform coding methods for 
data processing are realized in practica as transform of 
transmission, de-quantization/reconstruction. Another 
finite-impulse response digital ffiters where the mput 
rect calculation of convolution given by Fourier trans- 
,urier transform. For analog data, sampling and filtermg 
i speed of manipulation. 

provides a new method for sampUng analog or distal 
d to calculate new wavelet transforms as well as thwr 
JOS This immediately provides novel and efficient meth- 
ed on this new method of non^linear transform coding. 
veiet filters, the invention also provides new methods tor 
ansforms including the Fourier transform and its inverse. 
5vel aad efficient methods for digital filtering. By way 
tion, practical applications of the methods include non- 
,eech compression, speech recognition, speech synthesis, 
' for hearing aids, stiU two-dimensional image compres- 
m, video compression for purposes of telephony, precision 
Tonometry, denoising, interpolation, medical imaging, ge- 
of oU or other resources, and other emission/detection 
.tions, it may be necessary to store streams of input data 
sing or manipulation while concurrently coUectmg furthei 
r digital. It may also be useful in practice in paiticulai 
axcWve or otherwise store specific coefficients or constants 
>ry at run time rather than computing them at run time 

Description of Prior Art 

ihonormal bases of compactly supported wavelets. Com 
tlijTaihemaUcs 41 . pp. 909-996 (1988), Daubechie 

1 
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pre^l^ a class of wavelets which have been effectively ^PP^f-f.-^^'miiS^" 
SLts-. see. for instance. US Patents 5.453.945 by ^^'=^55,^,^:; ,f 
as w^ll J Daubechies' monograph 'Tten Lectures on VWets (J^f 
references. In another landmark paper, -An algorithm for the machine caicum, 
complex Fourier series", Mathematics of Computation 19, PP. 297-301 (1965), J-W. 

W J W Tulcey described a simplified algorithm for the calculation of Fourier 
Si^'te cTmp^- The utility of this fast Fourier transform FFT - w^^^^^ 
StrthoiSlds of important appUcations and i^Pj^^^^f f,^^^^^^^ 
pSents 3 871.577 by AveUax et al., 3,881.100 by Works et al., 4,051,357 by Borm^ot. 
I^rTheSBT. an owner's manual for the discrete Fourier ^^^^^"^ ' ^J^ii^Jg 
W Bri^s aid Van Emdea Henson and its references. A requiremaat for both the 
PPT me3 Sd fo7the multiscale filtering of the Daubecbies' wavelets is that the 
smpS^S iaTut data must be equally spaced though further methods of adaptive 
refinement of grids have also been developed in eadi context. 

Other examples of transforms whose efficient digital implementatiot^ ^^^P^^*^ 
SidTo which the invention is relevant include the ''^^f^'^^^^^t^. 
Bessel Laguerre, Hennite, Cbebyshev, Hotelling, Mersenne, and Fermat, see, lor m 
S U^Snts 3,891.443 by Lynch et al. and 4,093.994 ty Nussbaumer. 

Invention Summary 

Three new famiUes of wavelets, respectively called aHihm^Uc wavelets fan ^<^^^^^^^^ 
^^Z^ets are defined on the standard circle of radius one in the conaplex plane 
wtmpC^ue" function defined on this circle may be written ^^J^ially umqudy 
Cr combination of each family, and the approximation of specifi«^ mprnt 
data^TfiSte linear combination for eadi family leads to a correspond^g wavelet fi^^. 
E^imW^pSrupon a new quadrature on the circle. ^^^^^^ Th^Siput™ 
wSi relies ^ a novel method of non-equally-spaced samplmg. The o^*?^* 
^d^t filter plays the role in the Invention of an intermediate ^«P;^«»*^*^°^,f 
^Tdata and Ly one of a number of useful output transformations of data which 
may be computed from the intermediate quantity. 

With either analog or digital input data, it is convenient to think of the pv^ ito- 
Dut date as spatii data. The inverse wavelet transform or reconsirucUcm dgonthm 
r^uk^ c^culS spatial data from wavelet data and admits an especiaUy com.^ 
SrSip^entation: the values of the spatial function at certain ^X^/, 
l^l S^Lbitrarily fine spacing may be computed exactly using ^.JP^/^ ^ ^ 
ff wavelet coeffidents. This reconstruction algori^im «)mbineg with ^^J^^^' 
ter into a binary cascade giving an efiicient procedure for data compression wluch w 
2 ped^y SSlited to d£continuous input data and to iterative refi;^^^^* 
put dat^ JMrthermore, the methods extend direcUy to procedures J^^P'^^f 
of multi-dimensional data. Computer coding for the reconstruction algonthm itself 
?s sufficientiy abbreviated and low-level that It might be transmitted together with 
iSprS data, for instance, for down-loading from satelUte to home computer. . 

Jbr the calculation of classical transforms of spatial input data, it 
of the desired output data, for inatance, thelburier coeffici^ts. as frequency date. The 
transform from spatial to frequency date is accompUshed in ^-^^^^'^^^f^^^l 
of which approximates the input data by wavelets using the wavelet filter. The second 
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.tapis theo^lcuIationofFbuxierpodBdenuorother^^^^^^ 

Lid requires fcmiulM derived in previous *4 formules 

Brief Description of Figures 
PIG la illustrates the standard drde C of radius one in the 

-a 2)- 

FIGlblUustratesachordlnthetopsenJ-drdeofClebeledhyiiaalrixA- (° J). 

togetl^er vrtth two fartte "deseends.*" cb^ »« C ..beled ti.e matrb. p»d.oU 

J u rr r ^ ° 1 and T = ( M - The three chords 
C/7l and as indicated, where = I 1 I and ^ ij 

determine a triangle inscribed in C. an'd the Vertices of this triangle ax. indicated a. 
complex numbers. 

PIG ic illustrates a chord in'the bottom semi-circle of C labeled by a matrix A = 
^\ together with two further "descendaixt" chords of C labeled by the ma. 

trix products U^^A and T'^A as Indicated, where U'^ = i) ^ 

f 0 "i ) • '^'''^ '^^^'^^ ^ ^^^^ ^ ^' "''^ ""^^'"^ °^ 

this triangle are indicated as complex numbers. 

T7Tr 2 demcts a classical figure in mathematics called "the E^ey tesselation in Klein's 

^^dero?CbotT--t??. ^^^-^-n'b^'^'^rif 

by recursively taMng ^ies^^^ - -/^^ it^^S J brtS^I^pSg of 
ZuT^^Tc S^^td^in^ o^J^eS^^h: ^ey t^lation in increasing 
generation is a novel aspect of the invention disclosed here. 

PIG aa illustrates the circle C, the chord labeled I and the two triangles in the F^ey 

binations of cosine and sine functions, one such function next "*^Yp A 

dS^ed by the vertices, wh^e 9 is the usual a^igubr '^^'^^anate on t^ c^e. A 
fi!Sn ^7(60 i3 uniquely determined by the figure, where on -^ch such arc^ «c 
hm t^^lev^n% of the nearby combination of cosine and 
oiher parts Figures Sb-Se likewise determine analogous functions d^^) for other ma 
trices >1 =( " I). These functions taken together form the system of wavelets on 
which the invention is based, and they are muque to this invention. 
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FIG 3b illustrates the chord in the top semi-c^cle of C labded ^"-^lH^f^. ^^ 
Xe a > c. As before, the figure determines a function on the circle denoted 

FIG 8c Ulustrates the chord in the top semi-circle of C labded T^f^^^^^^J^ 
Xe .5 >- a. As before, the figure determines a function on the circle denoted iJx W- 

PIG 3d iUustrates the chord in the bottom semi-circle of C labeled ^^^^^^J^^^'f^.f: 
wh«e a > -c As before, the figure determines a ftincbon on the circle denoted iS^Cfl). 

PIG 3e iUustrates the chord in the bottom semi-circle of C }^'^^^^^7^''\^J^^^^ 
A where >= a. As before, the figure determines a function on the circle denoted 

PIG 4 illustrates the notation near the chord labeled by the matrix A = (^^ ^) » 
the top semi-circle of C; this notation wiU be useful in several subsequent discussions. 
FIG 5 provides a flow chart for the main driving routine in calculating the Fourier 
trEnsfomi' 

FIG 6 prOTddes a flow chart for the main recursive routine in calculating the Fburier 
transform. 

PIG 7 provides a flow chart for the procedure of updating coefiicients in calculating 

the Fourier transform. 

FIG 8 pro^ddes a flow chart for the procedure of processing terminal arrows of the 
recursion in calculating the Fourier transform. 

PIG 9 provides a flow chart for the main driving routine in calculating the inverse 
Fourier transform. 

PIG XO provides a flow chart for the main recursive routine in calculating the inverse 
Fourier transform. 

PIG IX illustrates the region U consisting of all complex numbers with pwitive imag- 
inary part. Map the disk I? to W via the function z^i{z + - 1); »f ^ ^ 
C^e the endpoints of a chord in Klein's model of the Farey t^elation as in Figure 2, 
then draw the semi-circle in U passing through the corresponding two pomte whi^^ 
perpendicular to the real axis. This produces another classical figure m mathematics, 
called the Tarey tesselation in the upper half-space model of hyperbohc geometiy" , 
which is indicated in Figure 11. 

Description of Preferred Embodiments 

In order to disclose the methods, it is necessary to develop some notation and termi- 
nology. Tb this end, consider the standard disk D of radius one about the wigm m 
the complex plane with its usual complex coordina.te z = x + %y, where » « v-i ^iia 
are real numbers, together with its boundary circle C aa illustrated in FiPjre la. 
An arrow is by definition an oriented chord of the circle C labeled by a two-by-two 
njatrix 5) , where the initial point of the chord is given by the complex number 
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Siia e C, and the final point of tbe chord is given by e C A special ^^<f c^^ 
doe (short for "distinguished oriented edge") is given by the diameter of C with 
endpoints -l.+l labeled by the matrix J J J) as in Figure la. Also required 

subsequently are the matrices 5=^5 "J)^) '^'^ = (n l)' ^ ^ (o l) ' 
any integer n, «hich satisfy the identities 5= = SUS - -T-\ STS = -U'^ of 
matrix products. 

Inductively define two separate binary cascades, one of top arithmetic arrows and an- 
other of bottom arithmetic arroiDs^ as follows: 

Basis: The doe, which is iUustrated b Figure la, is both a top arithmetic • 

arrow and a bottom arithmetic arrow. 

Top Induction: As in Figure lb, given the top arithmetic arrow labeled by 
yl := , there are also top arithmetic arrows labeled hy 

Bottom Induction: As in Figure Ic, given the bottom arithmetic arrow 
labeled by >1 = , there are also bottom arithmetic arrows labeled 

^ ^-'-=U. A)-^--'^=(°r V)- 

It is a classical theorem in mathematics that the chords underlying any two arithmetic 
arrows meet only at their endpoints, if at all, and that the family of all such chords 
decomposes the disk D into an infinite coUection of triangles with vertices m the circle 
C as is iUustrated m Figure 2; this figure is called "the Farey tesselation m Klem s 
model of hyperbolic geomet^y^ There is no suitable reference for this procedure m the 
literature, and it is for completeness that the expUdt recursive definition of arithmetic 
arrows has been given here. An arithmetic arrow is said to be of generation g if it anses 
from g appUcations of the inductive clause starting from the doe in case of either top 
or bottom arrows; the doe has generation zero by convention. 

Consider the infinite set Q £ C consisting of the endpoints of chords underlying all 
arithmetic arrows, and say a point of Q is of generation p if it is the endpoint of the 
chord underlying an arrow of generation g and no less; -1 and +1 are of generafaon 
zero by convention. There is a canonical enumeration of the generation g points of Q 
clockwise in C starting from -1 and the associated enumeration of Q itself by increasing 
generation. For example, the ordering begins with: 

-1^2f -2-t 2-i l-2i -l-3t 
-1, +1, +», -», zr+2i' ^T+i' 2 + i' l + 2i' -l + 3t' 

Given input data on the drde C. the sampling of it at the points QCC in this order 
of increaang generation will be referred to as the Farey quadrature. In practice, one 

5 
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might interpolate given data as required at the sampling points of the F^ey quadrature 
or restrict the F&rey quadrature to a circular arc in C. (In the subsequent discussion 
of mathematlcaJ basis, the saraple points of the F^ey quadrature will be seen to corre- 
spond to the rational points in the real line enumerated in increasing order, where the 
generation is the length of a corresponang continued fraction expansion.) 

Notice that the doe separates the triangle with vertices -1, +1, ~i from the triangle 
with vertices -1, +1, +i. In the same way, an arbitrary arithmetic arrow separates two 
triangles complementaiy to all the arithmetic arrows in !>• The endpoints of chords 
of these triangles in various cases are illustrated in Figures 3a-3e as functions of the 
matrix ( ° labeling the arithmetic arrow. Figure 3a depicts the doe itself 
with A = /, Figure 3b and 3c depict top arrows with o > c and c ^ a, respectively, and 
Figure 3d and 3e depict bottom arrows with a > -c and -c > a, respectively. The 
notation inside the circles indicates the endpoints of the chords as complex numbers, 
and the notation outdde the circles in each case of Figures 3a-3e is yet to be explained. 

Tb this end, define a Ixigonomeinc function to be a real-valued function 

/(fi) = C»cos5+ ^sind + 7. 

where $ is the usual angular coordinate on the unit circle C and o;, ^, 7 are real numbers. 
Next, define a piecewise trigonometric function on the circle to be a function f{Q) so 
that there is a finite collection of circular arcs J,- C C, for J - 1, .... J, any two of whidi 
meet at common endpoints in C, if at all, together with a coUection of trigonometric 
functions t^ifi) for j = 1, . . . , J, so that = tjie) if 5 e /,-. In other words, C is 
decomposed into a finite number of non-overlapping circular arcs, and on each such 
arc, / takes the values of a trigonometric function. 

There is a particular family of piecewise trigonometric functions, one such arithmeiic 
wavelet ^a(S) for each matrix A labeling an arithmetic arrow. These functions are 
defined in Figures 3ar3e, where the endpoints of the circular arcs are as described 
before, and each trigonometric function is written outside the circle next to its corre- 
sponding circular arc. There are thus five cases in the definition of arithmetic wavelet 
corresponding to the five cases in Figures 3ar3e ahready discussed, and in each case, the 
arithmetic wavelet itself is a piecewise trigonometric function defined using exactly four 
circular arcs. Arithmetic wavelets are a fundamental new ingredient of this invention 
and no doubt seem entirely ad hoc and unmotivated at this juncture of the discussion. 
The mathematical basis for them will be given later, but it is worth pointing out now 
that each arithmetic wavelet ^A(fi) is a once-continuously difierentiable function taking 
value zero at the points -1, +1, -i of the circle C, that is, for ^ « 0, -7r/2, -tt, except 
for i?t;-i (fl) and ^r~i (S) which take the value zero only at the points ±1. 

In order to explicate the definition of arithmetic wavelets, notice that any two-by-two 

matrbc ^ = ( c d) * self-mapping Ma iC-*C of the circle C as follows: 

/i-a _ (te4-&)-i(tc-K0 
^^[jTiTiJ - {ta + b) + iitc + d)' 

and in fact, the image of the endpoints, —1 and +1, of the doe under Ma are the 
endpoints of the arithmetic arrow with label A. Construct firom the function tf^Cfl) 
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defined on C in Figure 3a the corresponding vector field on C given by ^, where 
^ is the usual unit tangent vector field to C 

There is an explicit formula to he given presently for the transformation of suitable 
vector fields under the change of coordinates on C given by Ma^ Namely, the vector 
field = [or cos sin transforms under Ma to the vector field t^(fl)sgi 

where 



\ [2(c(i - 06)^ + (d' - + 7) + - - 7)] cos e 
(e) = | + [(od + l)c)/? + W(a + 7)-oc(tt-7)lsmfl j- 

+ I (2(cd + a6)j5 + (d^ + 6')(a + 7) - (a' + c')(o: - 7)] 

More generally, given a piecewise trigonometric function taking the values of the 
trigonometric function ij^Q) on the circulax arc 1} as before, the corresponding vector 
field /(fi)^ transforms under Ma to the vector field which on the circular arc MA\h) 
is given by t^W^. 

Changing coordinates on C by Ma using this formula transforms the vector field 
tfj(d) ^ to another vector field thereby defining the function 

Finally except for X = C^"^ and A = r-\ the corresponding arithmetic wavelet is 
given by 

■9a{0) = I'/lW - [o cos + /9 sin e+y], 
where ocJn are chosen so that tf(7r) = tf(0) - i5(-7r/2) = 0; the two special cases are 

du-i (6) = -du-i {$) + 2 sin e and ^t-« (^) = ''r-* (^) - 2 sin fl 

(which are defined in this manner to guarantee a symmetry of the upper and lower 
semi^drcles of C). This is the abstract definition of arithmetic wavelets which is made 
explicit in Figure 3. 

It is furthermore convenient to define 

so that for every label >1 on an arithmetic arrow, ?^(7r) = ?^(0) = i?A(-T/2) = 0. 
The arithmetic wavelets enjoy an important TtnormalizaUon property, namely, 

^BAiMAie)) = MMMs), 

where the prime denotes the usual derivative Mid where it is required that either both 
matrices B and A are matrix products of J7+Sr+i or both matrices S end A are 
matrix products of U~^, T~^. 

Just as a function / = may be expressed in its representation as a Eburier series / ~ 
'^ri cn^^^t SO too it may be expressed as a serie / w 5^ ^A^Ai^) of arithmetic wavelets, 
which together form a linearly independent set. The sum is over all arithmetic arrows, 
and the coefficients sa are called the ontAmettc wavelet coefficients. The calculation 
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of arithmetic wavelet coefficients from the input data / is called the arithmetic wavelet 
filter^ The arithmetic wavelet coefficients are not quite uniquely determined by / as 
will be explained. 

The reconstruction algoritlun or inverse wavelet transform provides a method for cal- 
culating function values at points of the circle from arithmetic wavelet coefficients. To 
determine the value of the sum ^Ca^a e-t a specified point q e Q Q C, draw a line 
segment from g to the origin in D\ this segment is contained in a finite collection of tri- 
angles complementary to all arithmetic arrows, whose boundary edges taken together 
form a finite set of chords; the corresponding finite set of arithmetic arrows A enu- 
merates the only coefficients ba which affect the value of 2^€S ^a^a at g- This is an 
important finiteness property of arithmetic wavelets which is crucial to the efficiency of 
the reconstruction algorithm. 

This completes the general discussion of arithmetic wavelets, the first family of new 
wavelets disclosed here, and the two further families of new wavelets are next defined. 

For each label >1 on an arithmetic arrow, there axe two fan wavelets 

n>0 

as well as two modular wavelets 

n>0 . 

It follows immediately from the definitions that there are the following relations among 
the wavelets for any label ^ on an arithmetic arrow. 

'^a(9) == <l>AW-4>UAi0) - M6)-4>USAi0) 

= ^£;-M(^)-2^^W+^a^(^) = ^U'^SAi6)-2ll>3AiO)+^USA(e). 

On the other hand as is not at all obvious from the definitions, (pi and are given by 
the following explicit formulas as piecewise trigonometric functions 



2 sin if 0 < < jt, 

-2cose - 2sind-2; if ir ^ 6 < 3Tr/2, 
2 cos 5 + 2 an e + 2; if 3ir/2 <9 <2n, 



and 



- 2 cos 5; If 0 < < TT, 
if < 5 < 27r. 



In fact, modular wavelets arise from V/ exactly the same manner that arithmetic 
wavelets arise from namely, i>j0)'§§ transforms under Ma to the .".Rector field 
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and-V'^Ce) = ^|>:^{e) - \a cos 6 + 0 sin 6 + y], where a,^,7 are chosen to 
guarantee that ^^(Q) = il>A{Tt) = ^>i(-7r/2) ^ 0; transforming ^j{6)'§§ using the 
change of coordinates given by Ms a likewise gives rise to ^SAi&) • The fan wavelets 
also arise from <f>j in the analogous manner. 

In contrast to arithmetic wavelets, neither the fan nor modular wavelets are linearly 
independent, and in each case, there is one linear relation for each arithmetic arrow, 
namely, 

<f>A — <f>UA = ^SA " ^'USAi 

On the other hand, the following collection of modular wavelets forms a basis 
{ipA ' A labels a top arithmetic arrow and c > a} 
U {tpsA i ^ labels a top arithemtic arrow and o > c} 
U {ipA ' ^ labels a bottom arithemtic arrow and - c > a} 
U {ipsA : -A labels a bottom arithemtic arrow and a > -c}; 

there is an analogous basis for the fan wavelets. 

Both fan and modular wavelets enjoy precisely the same finiteness property as arith- 
metic wavelets. Furthermore, both fan and modular wavelets enjoy renormalization 
properties, namely, 

^ba(Ma{B)) = M'M^BiS), 
where the requirements on -4,5 are as before. 

This completes the necessary abstract definitions- The methods will first be presented 
using arithmetic wavelets for real-valued input functions with extensions to modular 
and fan wavelets and complex-valued or multi-dimensional input/output data to be 
described later- As a point of notation throughout these discussions, there will be the 

. fa b\ 

tacit assumption that the entries in the matrix denoted A are given by ( . I - 



The Sburier Transform 

Included in the section entitled "Source Codes" is an implementation in the computer 
• language C of a preferred embodiment of the method described in this section employing 
arithmetic wavelets- 
Suppose that f{6) is some specified real-valued input fimction defined on the circle- 
The calculation of the Fburier transform of / proceeds by two main steps: 

Stepl : Choose a finite set S of arithmetic arrows and approximate / ^ 
Sx€S ^A^A calculating the coefficients Ca in a manner to be described in 
terms of the values of / at specified points using the flarey quadrature. 

Step 2^ Substitute into the expression in Step 1 the known Fburier expan- 
sions 19^(5) ~ e*^* in order to derive the approximation 

a^s 
9 
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The coefficients required In Step 2 are known theoretically and given for n ^4 o, ±1 by 



+ 2(c= + d=) 
-[(c + a)' + 



d- 
d + 

<^ + ^% + d) + zia + c)\ ■ 
Furthermore, the coeflScients c^, , c^x axe given by 



ir4i = eH{<?-d?± 2icd) + 9r[{h - df - (a - c)^ :f 2i(6 - d)(a ^ c)|/2 
+ d,(o2-62±2ia6) + ei[(6 + d)2-(a + c)='T2t(6 + d)(o + c)]/2, 

,rc^ = 2flfc(c3 + d=) - flr[(6'-d)* + (a-c)«] 
+ 25t(a» + &') - + df + (a + c)^], 

where the angles are 

. -1/ 2c6 . . ^ „w 2(& + d)(a + c) x 

(^^^ ^^-'^ ^(&-d)2_(a-c)2^* 

It may be useful in practice in particular applications to calculate and archive or oth- 
erwise store the constants required in the previous expressions for and recover them 
from memory at run tinae rather than computing them at run time. 

The Fourier coefficients co,c±i of / thus require special treatment which does not 
present additional challenges. For simplicify, suppose that the input function f{&) 
has these three coeflScients fixed lyy the condition that = 0 for = 0, —ir. 
More precisely, if f{d) is ^y real-valued function on the circle, then let /(fl) denote its 
norraalization defined by f \e) = f{B) + a cos 6 + j9 sin + 7, where Of, ^, 7 are chosen 
so that / takes value zero at 0, "'7r/2, The specified simplifying assumption on / 
thus amounts to the condition that / = /. 

The procedure for calculating wavelet coefficients is slightly different depending upon 
whether the corresponding arithmetic arrow is a top or bottom arrow. However, the 
antipodal map Q^it-^O transforms the bottom of the circle C to the top in such a 
way that it suffices to discuss only the case of top arrows. More precisely, the algorithm 
for bottom arrows arises &om that for top arrows by: 

• replacing the input data value f{d) with the input data value /(^ +'?r), 

• replacing the matrix ^4 in the described manipulations for top arrows 
by the conjugate matrix SAS for bottom arrows, and 

• replacing the trigonometric function ojcos 0 + /9sin 6 + 7 for top arrows 
by the trigonometric function -acos Q - j9sin (9 + 7 for bottom arrows. 

10 
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The algorithm and computer coding for bottom arrows is therefore derived without 
difficulty from that for top arrows, and the subsequent discussion is specialized without 
loss of generality to top arrows. 

In addition to the function f{d) defined on the circle, suppose that there is also specified 
a bandwidth of interest N (i.e. the method should calculate the Fourier coefficients Cn 
only in the range |n| ^ N) and a tolerance EPS (where contributions to the specified 
Fburier coefficients are regarded as negligible if they are of magnitude less than EPS). 
There are numerous variants of the method which depend upon other parameters which 
may be tailored to a given application, and some these variants will be discussed further. 

Steps 1 and 2 are merged into a single binary cascade which keeps track of an ongoing 
approximation to the Fourier coefficients as follows. A cascade element is called an 
arrow-structure and is defined to consist of the specification of 

• an arithmetic arrow A of generation ^, 

• the wavelet coefficient e^ on ^ the approximation of /, 

• the coefficients a, ^,7 of a trigonometric function defined as follows. The end- 
. points of the chord corresponding to A decompose C into two circular axes, 

actly one of which contains none of the points —1, +1, — i. The sum Y^b^a ^B^b 
ov&T all arithmetic arrows B of generation at most g except for A itself is given 
on this interval by the particular trigonometric function a cos + /? sin + 7. 

The coefficients a,j9,7 in an arrow-structure keep a lagged/updated running tally of 
the overall effect of what has come before it in the cascade; as a result of this technique, 
the method requires essentially no memory other than the storage of the running ap- 
proximation to Fourier coefficients and the stack required for the recursive computation 
in the cascade of arrow-structures. 

The next tedinical point involves a regulaiization in the calculation of arithmetic 
wavelet coefficients; this regularization is required beca\ise the approximation to / as a 
linear combination Y^j^^s ^a^a is not uniquely determined. Given an arrow-structure 
in the specified notation and referring to Figure 4, when sampling the input data at 
the point in the Earey enumeration with corresponding angular coordi- 

nate the one real input datum /(ffo) niust determine the two coefficients euA and 
eTA' According to the formulas presented in Figures Sa-Se, these new coefficients are 
constrained by 

^^^"^ - t(6 + d)3 + (a + c)=J +^ [(5 + d)2 + (a + c)«J +^ 
■ 4(ei/yt - sta) 8 Ba 3gn (a - c) 

'*'(6 + d)' + (a + c)2'^ (6 + d)2 + (a + c)2* 

where sgn (a — c) is a sign defined to be +1 if a > c and —1 otherwise. There remains, 
however, one real degree of freedom in the specification of euA, stAj and it is eliminated 
by demanding the best least-squares fit to the values of / at the next generation of points 
iaKS-'ia^t°i ' iS^-lSSi with respective angular coordinates ^i.fla, as illustrated 
in Figure 4. Updating the lagged trigonometric coefficients in accordance with the 
formulas presented In Figures 38p3e produces trigonometric functions Vi(j9) and VziS) 
which give the values of the ongoing approximations at ^i.and 02, respectively, and 

11 
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whose coefficients are explicit inliomogeneous linear functions of evAt^TA, namely, 

Vf(8i) - vie + viu evA + v»t crAi for t = 1, 2. 
The coefficients eie expUcitly described most concisely in pseudo-code, where e = e^, 
alp = a, bet = /3, gem = 7. bpd = ft + d and ape 1= a + c: 

if(a > c){ 

axl=alp+(d*d-c*c+2*(a*c-b*d)+a*arb*b)*2*e; 

ayl=2*(b+b-a*a); 
bxl=bet+(c*d-a*d-b*c-a*b)*4*e; 

byl=4*a*b; 

cxl=gam+(2*(a*c+b*d)-c*o-d*d+a*a+b'*b)*2*e; 

cyl=-2*(a»a+b*b); 

ax2=alp-|-4*(d*d-c*c)*e; 

ay2=2*(c*c-d*d); 

boc2=bet+8*e*c*d; 

by2^4*c*d; 

cx2=6am-4*e* (c*c+d*d) ; 
cy2=2*(c*c+d*d); 

> 

axl=alp+4*e*(a*a^b*b); 
ayl=2*(b*b.a*a); 
bKl=bet-8*e*a*b; 
byl=4*a*b; 

cxl=!gam-|-4*e*(a*a+b*b); 
cyl=-2*(a*a+b*b); 

a3c2=alp+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)»2*e; 
ay2=2*(c*c-d*d); 

bx2=bet+(a*d-a*b+b*c+c*d)*4*e; 
by2=-4*c*d; 

cx2=gam+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*a; 

cy2=2*(c*c+d*d); 

} 

vlc=cxl+(((b+bpd)*(b+bpd)-(arfapc)*(a+apc))*axl 

-|-2*(b+bpd)*(aH-apc)*bxl)/((b+bpd)*(b+bpd)+(a+apc)*(a-|-apc)); 

vlt=cyl+(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*ayl 

-|-2*byl*(b+bpd)*(a+apc))/((b+bpd)*(b+bpd)-l-(a+apc)*(a+apc)); 

vlu=8./((b-l-bpd)*(b+bpd)+(a+apc)*(a+apc))j 

v2c=cx2+(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ax2 

+2*(d+bpd)*(c+Rpc)*bx2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

v2u=cy2+(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ay2 

+2*(d+bpd)*(c-i-apc)*by2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

v2t=-87((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

It is easy to minimize [/(^i) - Vi(fli)]' + [/(fa) - ^3(^2)]'' 93 a function of euA^tTA 
subject to the given constraint, and this procedure provides a suitable regulanzation. 

12 



PAGEIM'RCVDAT7/6I2I)089:37:3I AM [Eastern DayfigtitTiine]'SVR:USPTO{Fn 



08-Jul-06 07:36311 Froir 4142715770 T-S20 P. 014/074 F 



Weighted least-squares fits or other statistical treatments using input data values at 
further generations allow alternative regularizations. 

Two other simple regularizations which have also been useful are: to require that 
er;^ - 0 if c > fl and ejA ^ 0 if o > c, or, to require that euA/sTA = 
d)]/[c(Q+c) + d(6+d)l; the latter regularlzation is related to imposing differentiabUxty 
of the approodmation at 9q. 

The final technical point involves the stopping criteria and terminal processing for the 
cascade of arrow-structures. The basic stopping parameter is NVAN, where a branch 
of the cascade terminates whenever there have been NVAN consecutive generations 
of cfcpring whose contributions to all coefficients in the bandwidth N have been 
of magnitude at most EPS. There are further parameters governing stoppmg cnteria 
which may be introduced depending upon the scale and resolution of features in the 
input data, including a minimum or maximum number of generations, the size of the 
angle subtended by the chord of a terminal axrow-structure, and so on. 

Another useful stopping criterion is to require that |ex|J|^>i|l be negligible for sev- 
eral generations, where W^aW denotes the i^.norm of i?x for which there is the a 
prion estimate < I0|«c + bd]-^/*; that is, one terminates the cascade when 

NVAN consecutive generations of coefficients Ba have satisfied the condition that 

When NVAN generations of offispting contribute negligibly to the specified bandwidth, 
the cascade is terminated with those offepring of greatest generation. Suppose that A 
labels the arrow of such a terminal offspring, adopt the notation in Figure 4, and let 
JQC denote the circular arc with endpoints ^ which contains none of the 
points -1, -i, +1 6 C, The final update alpi cos * + betf sin ^ + garni, for t = 
1,2, of the running trigonometric function on J in accordance with Figures 3ar3e is 
prescribed in pseudo-code with notation as before as follows: 

if(a > c){ 

alpl=alp+e*2*(a*(a+c)-i-d*(d-b)-b*(b+d)-c*(c-a)); 

betl==bet+e*4*(c*(d-b)-a*(d+b)); 
gaml=gam+e*2*(a*(a+c)+c*(a-c)+b*(b+d)-d*(d-b)); 

alp2=alp+e*4*(d-c)*(d+c); 

bet2=bet+e*8*c*d; 

gam2asgam-e*4*(o*c+d*d); 

} 

else{ 

alpl-alp+e*4*(a^b)*(a+b); 

betl=bet-e*8*a*b; 

gaml— gam+e*4*(a*a+b*b); 

alp2=alp+e*2*(a*(arc)+b*(d-b)-c*(c+a)+d*(b-|-d)); 

bet2=bet-|-e*4*(a*(d-b)+c*(b+d)); 
gam2=gam+e*2*(a*(a.c)-c*(a-Fc)-b*(d-b)-d*(b+d)); 

> 

Now serially letting t denote the trigonometric functions n and ra, again adopt the 
notation of Figure 4, where A labels the arrow corresponding to t. The updated r 
is compared with another trigonometric function a defined by extrapolation which is 

13 
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uniquely determined by demanding that it agree with / at the three angles comprising 
the next two generations, 6 = 60,0^,62 in the notation of Figure 4. The Fbun^ 
coefficients of the piecewise trigonometric function which takes values '^(«) "/Cf) on 
the circular axe J and value zero otherwise are easily computed using standard formulas 
and added to those of the ongoing approximation. 

Other algebraic or statistical samplings of higher generation offspring allow alternative 
final processings. 

Fbr histanoe, one xnight save computational expense and simply terminate the cascade 
with no final processing at all. 

Furthermore, the terminal arrow-structures could be stored for restart or iterative 
refinement capabilities. 

One might alternatively record in an arrow-structure different transforms of the trigono- 
metric function defined by a,/3,7. For instance, it is useful to take advantage of the 
renormallzation property of wavelets and change these coefficients by trar^orming the 
vector field [a cos B+p sin 0+f]-§s by the change of coordinates This tecluuque 
avoids the necessity of calculating certam large integers attendant to the calculation of 
wavelet coefficients as will be illustrated later. 

One favorable aspect of the method is its advantageous mix of floating-point and integer 
operation types. Moreover, except for the coefficients co, c±i, the implementation ofthe 
algorithm is purely algebraic, that is, requires only addition and multipUcation. Rur- 
thermore, by its very nature as a cascade, the algorithm is amenable to paralleUzation 
and efficient hardware implementation. 

This completes the description of the implementation of the method of calculating 
Fourier coefficients using axithmetic wavelets. 

Several points will next be addressed which are special to the implementation of the 
methods using modular wavelets. Again there is a basic Step 1, the c^culation of 
the modular wavelet transform where the input function is approxmiated as 
T:^gAi>Ai9), and a basic Step 2, substitution of known expressions for the Fourier 
^efficients of the modular wavelets; these two steps are again conveniently merged 
into a binary cascade of arrow-structures in practice. 

To elucidate Step 1 for modular wavelets, refer again to Figure 4 for the notation near 
the axithmetic arrow labeled by the matrbc A. The single value of the input func- 
tion fieo) together with the ongoing calculation ofthe updated trigonometric function 
o cos e sin e +7 in the arrow-structure this time uniquely determines the modular 
wavelet coefficient 

_ /(gp) - cos go + j3 sin gp + 7) 

since ^Bida) = 0 for every labeling matrix B of generation greater than A] there is 
thus no need for regulariaation with modular wavelets. 

As to Step 2, recall that a basis of modular wavelets was given before by a collection 
V;fl(g) of wavelets where A is the label on an arithmetic arrow and either 5 = A or 
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B = SA. In each case, one can calculate that ^si^) has Jbnrier expansion 



An9 



|rt|>i 



where 



2 (o2 + 6^) 



VS+ic) Kb+iaJ 



; i£B = SA, 



\s4 



and 



27rdS 



27rdf 1 = 2ndf 



2(a2 + 6»)[tan-i3l^-tan-i3|^]; i£B = SA, 



Thus, the method of calculating Pburier transforms using modular wavelets is easily 
derived from the method using arithmetic wavelets. In feet, it is a simple matter to in- 
corporate the two minor modifications just described and alter the included source code 
employing arithmetic wavelets to produce source code employing modular wavelets. 

In the same way for fan wavelets, there is again a Step 1 and Step 2, which are merged 
into a binary cascade of arrow-structures. As for arithmetic wavelets, Step 1 for fan 
wavelets again requires a regularization scheme. As to Step 2, the formulas for Fourier 
coefficients of fan wavelets can be derived from those given before for modular wavelets 
using the identity <pAi&) = - ^ua{B) mentioned before. Again, the included 

source code is easily modified to employ fan wavelets. 



Flow Charts for the Fourier 'transform 

Flow charts for various procedures used in the calculation of Fourier transforms are 
presented in Figures 6-8. .For definiteness, the flow charts and source code refer to 
arithmetic wavelets without renormalization, but the renormalizing source code will 
also be included for completeness. For purposes of brevity in these flow charts, the 
arrow-structures defined before are lefetred to in these figures dmply as "arrows". 

In Figure 5 is presented a flow chart for the m^ driving routine. The input data is 
normalized as indicated in program segment 1. There are separate recursions estab- 
lished in program segments 2 and 3 for the top and bottom of the circle respectively. 
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Each recursion U established with a call to the subroutine gener, which then 
ta^n. The Fourier coefficients are stored internally with an overall suppressed factor 
S TT and output data normalization of multiplying by tt is accompUshed m program 
segment 4. Output Fourier coefficients are finaUy displayed before eating m progr^ 
5 Pro^ segments 2 and 3 are entirely independent and could be performed 
iu^ld; mor^^nerSy, ea^Ji of program segments 2 and 3 could be decomposed 
further into multiple parallel procedures. 

A flow chart for the main recursive routine gener is given in Figure 6. The procedure 
starts with a test in program segment 6 to determine ifc 

• the argument envy >= NVAN, in which case there have been NVAN 
consecutive generations of negligible contributions, and 

• the argument generation >= MING, in which case there have been a 
required minimum number of iterations in the recursion. 

If both of these inequalities hold, then the procedure passes to progr^ se^ent 15 
where the cascade is terminated and the output data corrected with a call to the 
subroutine prune. 

In the contrary case that the recursion should continue, the procedure passes to prosram 
segment 7, where the descendant arrows in the cascade are determined usmg the lea^ 
8q?^es fit to the next generation of data as described before to compute the regu arized 
wJtdet coefficients of the descendants. These are combined with integer calculations 
to update the lagged trigonometric functions. The procedure then passes to pro^ 
segment 8, where the ongoing appnsdmations to Iburier coefficients are updated to 
in^dc coitributions from the two descendant arrows with a call of the subroutine 
Charles for each descendant. The procedure continues with a test in program segment 
9 If the contribution calculated in the subroutine charles for either descendant to any 
Eburier coefficient in the bandwidth was non-negUgible, then the recurtive argument 
envy is set to zero in program segment 10. In the contrary case that both oontnbutions 
to i Fourier coefficients in the bandwidth were negligible, the control parameter envy 
is increased by one In program segment 13. 

In any case, the procedure passes to program segment 11, where there is a test on the 
number of generations. In case too many iterations of the recursion have occurred 
it is terminated in program segment 12. and suitable warning of non-convergence of 
the method is written to the output. In the contrary case that the maxmium number 
of generations has not been reached, program segment 14 establishes the recursion by 
calling gener once for each of the descendant arrows with the updated control parameter 
envy and an incremented generation. 

In Figure 7 is presented a flow chart for the subroutine charles. Calling the subroutine 
Charles has the effect of updating the ongoing approximations to the Iburier co^aents 
for a single argument arrow and returning a flag which keeps track of whether a^ 
such contribution has been non-negUgible. The flag is cleared in program segment 16, 
and several preUminary calculations are accomplished in program segment 17. ifte 
OtOCedure passes to program segment 18, where it ia determined whethM to calci^te 
the 0,+l,-l Fourier coefficients. If these are to be calculated, then the procedure 
passed to program segment 19, which calls a subroutine to update these three Fourier 
coefficients. Program segments 20 and 21 set the flag as required dependmg upon 
whether these three contributions are non-negli^ble. 
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In any case, the procedure passes to program segment 22. If the loop ^^e^^^^^^ 
width is complete, then aU Fourier coefficients have been updated so the subrDutine 
Charles returns the flag in program segment 23. In the contrary case, ,^7^* 
coefficient is updated in program segment 24. The flag is set as ^^^^^^^^ ^depend^ng on 
whether the current contribution was non-negligible in program segments 25 aad 26 
The iteration over bandwidth is established by the update m program segment 27 and 
the return to program segment 22. 

In Figure 8 is presented a flow chart for the subroutine prune, which modifies the u- 
ray c5 Fourier coefficients when terminating the cascade. The final update oflaggwi 
trigonometric functions to produce n, is performed in program segment 28. "^ere ^ 
one such function for each possible descendant of the argument arrow. The procedure 
passes to program segment 29, where further input data is samp ed m order to com- 
pute two ^gonometric extrapolations cri.aj, one such extrapolation for each po^ible 
descendant of the argument arrow. Each function - , for a = 1, 2, is ^^<^^ ^ 
described before, and the Fourier coefficients of the truncated functions are «lded to 
the ongoing approximations of Fourier coefficients in program segment 30. f^^^ 
segments 31 and 32 implement the calculation of 0,+l,-l Fourier coefficients if desired, 
and the subroutine prune is terminated with the return in program segment 33. 

The Inverse Fourier Transform 

Included in the section entitled "Source Codes" is an implementation in the computer 
language C of a preferred embodiment of the method described in this section employmg 
arithmetic wavelets. 

Suice the method for calculating the inverse Rjurier transform is similar in ^irit and 
execution to that for the Fourier transform described in detail before, it need only be 
briefly discussed. The input data includes the specification of a collection of Pouner 
coefficients c„ in a given bandwidth N, 

There are again two main steps in the calculation: 

Stepl : Calculate arithmetic wavelet coefficients ca from the given Fourier 
coefficients c„. 

Step 2: Use the wavdet coefficients from Step I and the reconstruction al- 
gorithm to output values of the function ca^aW at the Farey quadra- 
ture points on the circle. 

For Step 1, there Is again an explicit formula firom previous theoretical work. Deep 
mathematical results prove that for nj^O, ±1, 

e**^ aa cos + 1 sin 7i5 

where the sum is over all arithmetic arrows, the underlying chord of which has endpoints 
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^,7/ e C, and 



-1, ns 0(4)1 

0, nBl(4)i 

-1, ns2(4); 

0, ns3(4)i 



Cj — 





nsO(4); 




nsl(4); 




n = 2(4); 


I 0, 


n = 3(4); 



I*" 



0, n 5 0(4) 

0, n = 1(4) 

-i, n = 2(4) 

-1, n3 3(4) 



Substituting this expression for e*"^ Into the given Fourier expansion f($) ~ J^n 
gives 

A 



where 



and 



This expression /(e) ^d^ + cie<^ + cLie"'^ + ^aC^) of /(fl) constitutes Step 1 

in the calculation of inverse Fourier transform. 

In fact, the calculation of wavelet coefl5cients in Step 1 can be implemented using 
mostly integer operations to improve run-time dramatically. To see this, notice that 
the formulas ^ = = 3^ for the endpoints of the chord underlying the arrow 

labeled by the matrix = d) substituted into the expression for ca to 

yield 



where 



+ ac + i)"" + {bd + ac - i)"} + + ac){ (6d + ac + i)"" - (W + ac - »)"} 



In the obvious implementation of this formula to calculate using integer operations, 
the intermediary integer expressions grow too large to be practical for n large if the 
arrow corresponding to A is of large generation. 

This difficulty is overcome by wavelet renormalization, where these attendant large 
integer calculations are obviated entirely as is illustrated in the included source code. 

Step 2 in the method for calculating inverse Fourier transforms depends upon the 
reconstruction algorithm to output function values on the circle taken by the Imear 
combination f(9) « cj, + c^e*^ + cLxfi"** + Za^a ^Ai9). The finiteness property of 
the reconstruction algorithm allows the exact calculation of these function values at 
the sample points Q C C in their ordering determined by the Farey quadrature. The 
two steps are again conveniently merged iato a single binary cascade as follows. 
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There i3 only one control parameter, which is called SCALE, where the method is 
required to determine at least one output value in each circular arc subtending an 
angle SCALE. Thus, 85 SCALE is decreased the grid of output values on the circle is 
refined. A binary cascade is established using arrow-structures which is similar to that 
for the fburier transform before. However, the wavelet coefficients are given in this case 
in closed form by the formula for e^A (so no least-squares fit is required). Because of the 
finiteness property of the reconstruction algorithm, exact output values are determined 
at each stage of the recursion. The cascade is terminated with an arrow whenever the 
chords underlying both descendant arrows subtend angles less than SCALE. Output 
data can be written to an array either when calculated to store data in the ordering of 
the Farey quadrature or when terminating the cascade to store data in the standard 
ordering on the circle. Furthermore, the terminal arrow-structures could be stored for 
restart or iterative refinement capabilities. 

This completes the description of the implementation of the method of calculating the 
inverse Fourier transform using arithmetic wavelets. 

Several points will next be addressed which allow the extension of the methods to 
modular and fan wavelets. Again there is a basic Step 1, the calculation of wavelet 
coefficients from Fourier coefficients, and a Step 2, the reconstruction algorithm; in 
each case of fan or modular wavelets, these two steps are again conveniently merged 
into a binary cascade of arrow-structures in practice. The finiteness condition on fan 
or modular wavelets mentioned before renders Step 2 entirely an^ogous to that for 
arithmetic wavelets. As to Step 1, simply substitute the identities '^a-^a- ^va = 
^ijA - SV'i^ -h i>U''^A described before into the expression above for e^"^ in terms of 
arithmetic wavelets to immediately derive corresponding expressions in terms of fan 
and modular wavelets. 

Thus, the method of calculating inverse Fourier transforms using fan or modular 
wavelets is easily derived from the method using arithmetic wavelets. In fact, it is 
a simple matter to incorporate the minor modifications just described and alter the 
included source code employing arithmetic wavelets to produce source code employing 
fan or modular wavelets. 



Flow Charts for the Inverse Fourier Transform 

Flow charts for various procedures used in the calculation of inverse Fourier transforms 
are presented in Figures 9-10. For definiteness, the flow charts and source code refer to 
arithmetic wavelets without renormalization, but the re-normalizing source code will 
also be included for completeness. For purposes of brevity in these flow charts, the 
arrow-structures defined before are referred to in these figures simply as "arrows" . 

In Figure 9 is given a flow chart for the main driving routine. Modified Fourier coef- 
ficients Co,ci,cli are computed in program segment 34. The formula given in Step 1 
for the inverse Fourier transform is applied in program segment 35 to calculate the 
arithmetic wavelet coefficients in terms of Fourier coefficients for a particular family 
of nine arrows. These nine arithmetic wavelet coefficients are required to initialize the 
trigonometric functions for recursions established in program segment 36 with calls 
to the recursive subroutine gener. Recursions with initial conditions corresponding to 
^ = £/,T,r-^, C/'"^,r"Sr"^C/-^, are undertaken in program segment 36. These 
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recursions, as well as further possible subrecursions derived from them, are indepen- 
dent and may be undertaicen in parallel. This elaborate driving routine is convenient 
because e^^^ is expressed in Step 1 ofJ;he inverse Fourier transform method in terms 
of the normalized arithmetic wavelets which differ from the arithmetic wavelets i9a 
themselves only for A = £/"^, T^^. 

The procedure passes to program segment 37 to perform an overall correction of the 
function values usiog the modified Fourier coefficients co.<i<^-i- Output data is dis- 
played and the procedure terminated In program segment 38. 

In Figure 10 is given a flow chart for the main recursive routine gener. In the cas- 
cade, the two descendant arrows of the argument arrow are generated in program 
segment 39 employing Step 1 of the inverse Fourier transform method to calculate the 
wavelet coefficients of the descendants. These expressions are then used to update the 
lagged trigonometric functions in program segment 39. A test is performed in program 
segment 40 to determine if the chord underlying each descendant arrow subtends a 
sufficiently small angle, as estimated by the comparison of an integer quantity with 
SCAT; the relation with the control parameter SCALE discussed before Is given by 
scat' = [exp sinh"*(l /SCALE)] ^ In case one of the chords subtends too large an 
angle, the procedure passes to program segment 41. The recursion is estabHshed in 
program segment 41 with two calls to gener, where the arguments are given by the 
two current descendant arrows. In the contrary case, the procedure performs the final 
update of the lagged trigonometric functions in program segment 42, then stuffs the 
output array with the new function values in program segment 43, and finally returns 
in program segment 44. 



Data Compression 

As with any method for data compression based on transform coding, five basic ma- 
nipulation are required, namely: 

• wavelet filtering, 

• quantization, 

• storage and retrieval or transmission 

• de-quantization, and 

• reconstruction. 

The wavelet filter has aheady been fully disclosed as Step 1 for the calculation of the 
Fourier transform, and the wavelet inverse filter or reconstruction algorithm has likewise 
aheady been described as Step 2 for the calculation of the inverse Fourier transform. 
Application-specific quantization is done according to psychovisual or psychoacoustic 
thresholds. Quantization and storage are furthermore merged with wavelet filtering 
into a single binary cascade as before, and retrieval or transmission are likewise merged 
with reconstruction into another single binary cascade. The reconstruction algorithm 
is exact. Furthermore, source code for it may be transmitted along with compressed 
data since the coding is sufficiently abbreviated and low-level. 

As is well-known, useful compression by transform coding requires an efficient specifi- 
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cation of the basis elements. This is especially advantageo\is for the new wavelets: a 
top arrow of generation g labeled by the two-by-two integral roatriitf A is specified by 
0 bits, namely, by writing A uniquely as the matrix product of factors U or T. For 

(13 9\ 
^ j factors as 

(lo 7)'={o 00 l)(o l)0 0' 

SO the coefiScient for the wavelet 4 is coded by the binary sequence TUUUTTU . 

The method is also espedaJly well-suited to progressive pictxure build-up or other it- 
erative refinement as follows. The Farey quadrature determines a linear ordering 
on the set Q of quadrature points in the circle C as before. The ordered subset 
Q\{-i,±l} C Q C C is put into bijective correspondence with the collection of all 
arithmetic arrows (where \ denote the set-theoretic relative complement); as illustrated 

in Figures la-lc, the arithmetic arrow labeled by -A ^ ^ ^ corresponds to the com- 
plex number 

^g±Ms+g6C, foratoparrow; 

\ j^ljjtjjolej € for a bottom arrow; 
^ t e C, for the doe. 

Thus, the ordering from the Earey quadrature on ±1} determines a conespond- 

ing ordering on the set of all arithmetic wavelets. Storage, retrieval, transmission, and 
reconstruction of wavelet cpeflBcients is accomplished at a specified input or output data 
resolution in this canonical ordering since energy compacts to the lesser wavelet coef- 
ficients. In the case of real-time compression and transmission or other manipulation 
not requiring retrieval, the method can be further improved by the manipulatioa of 
previously computed wavelet coefficients in parallel with the calculation of subsequent 
ones. 



Complex and Multi-Dimensional Data 

There is nothing remarkable about the extension of the invention to complex-valued 
input data = + % where u and u are real-valued functions on the circle 
C7. Concentrating for definiteness on arithmetic wavelets, one simple approach is to 
separately and independently transform u and t;, 

u(e) ^Y^tA and ^Y.^a ^a{B\ 

using the techniques already described in order to calculate 

Another direct approach is to allow complex values for the input function / and solve 
for complex wavelet coefficients using the same algebraic formulas as before. 
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Turning now to multi-dimensional input/output data, suppose that there are JZ > 1 
input functions each of which has M > 1 independent variables, and define the torus 
to be the M-fold Cartesian product 

M times 

Coordinates on are given by M-tuples 5^= ^ji • • • i ^Af ) of angles, and the input 
data is specified by an iZ-tuple of multivariable real- or complex-valued functions 

F{9)^U\S),ne),.,.j\e)). 

Using the buection described in the previous section on data compression, each quadra- 
twe point Z 6 uniquely determines a corresponding arithmetic wavelet 

^z{0)\ by convention, one also defines ^~i{6) = = 1. Concentrating on arith- 

metic wavelets for definiteness, define the corresponding mulUwavelet 

As before, the input function F must be normalized- To this end, define a trigonometric 
monomial to be a product of the form X^{ei)X^{e^) - • ■ X^^Om), where each X'{9) is 
given by one of = cos I?, = sin fl, or X'(ff) ^ 1. The values of F at 

the points in {— i, ±1}^ uniquely determine coefficients in a linear combination u0) 
of trigonometric monomials so that = F(ff) for each point of {-i, ±1}'^. Just 
as in the one-dimensional wavelet filter already discussed, now one must consider the 
normalization F{i) = F(d) - u{§), so is zero at each point of ±1}^- 

Adopting the Farey quadrature in each factor circle C of the torus C^, there is an 
induced lexicographic ordering on the set £3^\{-i,±l}^ C of potential sample 
points Z 6 C^. One separately approximates each coordinate function 

S 4 ^Zi^^ for i - 1, 2, . . . , i2, 

of ^(d) as a finite linear combination of multiwavelets with application-specific zonal 
sampling for optimum energy compaction; this calculation may employ least-squares or 
other regularization in each factor as well as renormalization as for the one-dimensional 
wavelet filter. The methods are again elegantly implemented in practice as M nested 
binary cascades. The final processing for the cascades of arrow-structures in the cal- 
culation of multi-dimensional Fourier transforms may involve adding suitable linear 
combinations of trigonometric monomials based on further sampling, or one may sim- 
ply truncate with no final processing at all, as in the one-dimensional case. 

These same remarks apply verbatim to fan and modular wavelets. 



Mathematical Basis 

**The decorated Teichmuller space of punctured surfaces", Communications in Mathe- 
matical Physics 13, pp. 299-339 (1987), written by the author of this patent applica- 
tion, began the study of certain abstract geometric coordinates called lambda lengths 
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in the context of hyperbolic geometry on Riemann surfaces. The ideas developed in 
this and other related papers have found applications in high energy physio. A more 
recent publication in this series of about 20 papers is entitled "Universal constructions 
in IfeichmuUer Theory", Advances in Mathematics 98, pp. 14a-2lS (1993). Here it 
was shown how to generalize lambda lengths to provide coordinates on a certc^n space 
of homeomorphisms (i.e. continuous byections whose inverses are also continuous) of 
the circle, where there is one generalized lambda length for each arithmetic arrow. 
In further recent work described in '*The Lie algebra of homeomorphisms of the cir- 
cle", Advances in Mathematics 140, pp. 282-322 (1998), written jointly with Feodor 
Malikovj there is described a repr^entation of the coordinate deformations of these 
generalized lambda lengths as explicit vector fields on the circle. These vector fields 
axe called "normalized elementary vector fields", and their study led to other vector 
fields on the circle called **fans" and "hyperfans". These three families of vector fields 
correspond, respectively, to the arithmetic, fan, and modular wavelets described here 
via the identification of the vector field /(5)^ with the function f{0)\ there was, how- 
ever, no discussion of wavelets, sampling of data, approximation of functions, or any 
algorithms in the published literature. 

Each arithmetic wavelet is once-continuously diflferentiable on the circle, compactly 
supported, and locali2ed in space. Arithmetic wavelets are not compactly supported 
in frequency, but instead, the frequency profile of an arithmetic wavelet is given alge- 
braically by the formula for used in Step 2 in the calculation of Fourier transforms; 
a non-compactly supported localization in frequency follows directly from this. The 
asserted formula for can be derived without much difficulty but in several cases 
directly from Figures 3ar3e integrating twice by parts the usual expression for Fourier 
coefficients using the fact is once-continuously differentiable. More difficult to derive 
is the formula given for the exponential functions e*"^ in tenns of arithmetic wavelets 
which was used in Step 1 of the calculation of inverse Fourier transforms; this subtle 
computation with lambda lengths is given in 'The Lie algebra of homeomorphisms of 
the circle" . 

Similar remarks apply to fan and modular wavelets, but fan wavelets are only contin- 
uoixs (they are not differentiable), and modular wavelets are not even continuous. 

In order to explain the failure of orthonormality of arithmetic wavelets and illuminate 
the regularization used in the arithmetic wavelet filter, fix some matrbc A labeling 
an arithmetic arrow whose imderlying chord has initial point g in the circle C, The 
sequence U^A of matrices label the arithmetic arrows whose underlying chords also 
have q as initial point, where n denotes an arbitrary integer. There is one infinite 
linear dependence 

n 

for each q e Q, and this explains the additive degree of freedom in the calculation of 
ei/A and erA m the wavelet filter which is handled by regularization as was discussed 
before. This additive degree of freedom is also manifest in the initial specification 
e/ 0 of wavelet coefficient on the doe in the wavelet filter. 

In order to better understand the Parey quadrature, it is convenient to transform the 
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disk D to the upper half-plane W = {tt + iv : v > 0} by means of the transform 

. z + 1 
z — \ 

This function K \3 a. homeomorphism on the interior of D which maps the unit circle 
C bijectively to the real-axis {u + iv:v = 0} together with an additional point oo 
at infinity, where Jir(-l) = 0, K{-i) = 1, and = oo. Rirthermore as g € Q 

varies over all quadrature points, the real number K{q) varies over all rational numbers, 
where the rational number E corresponds to the quadrature point 2=^ e <5. In order 
to Illustrate the Farey tesselation itself in W, whenever 91,92 S Q are the endpomts 
of the chord underlying an arithmetic arrow in £>, draw the semi-circle in U with real 
endpolnts 91,3-.; in the degenerate case that one of the quadrature points is infinite, 
say the point 92 = 00, then draw the vertical ray in U with endpoint q\. Applying this 
procedure to each arithmetic arrows produces the classical model of the Farey tesselation 
in U that is indicated in Figure 11. This figure iUustrates the relationship between 
the IWey tesselation and the indicated circle packing, where the circle in the figure 
that is tangent to the real-axis at the rational point p/q has diameter 9- . For further 
information on number-theoretic aspects of this classical figure, see, for instance. Ford's 
book "Automorphic Functions", Chelsea (1972). 



Source Codes 

In order to adequately disclose the methods, it is the purpose of this paragraph to 
include source codes for their implementations in the computer language C. Source 
code is presented for each of the following two Implementations. 

• awft is the implementation of a preferred embodiment of the method 
for computing the Fourier transform of real-valued input data defined on 
the circle. The algorithm depends upon a bandwidth N, tolerance EPS, 
generation cut-off NVAN, and minimum generation MING as described 
before. A further flag SM0DB~1 determines that the output data will 
include the 0,+l,-l Fourier coefficients, and SMODE=0 otherwise. (The 
only configuration in either code which requires any library routines is 
SMODEs=l, which requires math.h.) 

• awift is the implementation of a preferred embodiment of the method 
for computing the inverse Fourier transform of specified complex-valued 
Fourier coefficients defined in a specified bandwidth N. The output data 
is controlled by a single parameter SCAT, which is related to the param- 
eter SCALE discussed before by SCAT=[exp sinh " (-1)(1/SCALE)] * 2. 
(SCAT varies inversely with SCALE and may be compared with conve- 
niently derived approximate quantities.) 

The comparison of the subroutines tgener in awfc and in awift illustrates the easy 
transition between source codes for real and for complex data. The corresponding 
complex-awft and real-awift are therefore easily derived from the included source codes. 

24 



PAGE 2SI74'RCVDAT7n69:37:31 AM [Eastern DayligtitTiine]'8m-USi^^ 



06-Jul-06 07:4030 From- 4142715770 T-520 P.0Z6/074 F- 



Furtherroore, the method for compression of one-dimensional input data amounts to 
the first step of awft, then quanti^ation/de-quantization, then the second step of awift. 
The included source codes together thus also implement the method for compresaon 
of one-dimensional data since source code for it may be extracted from the included 
source codes. FinaUy. it is straight-forwaxd to write the further driving routine required 
for multi-dimensional data compression, so the source codes included with this patent 
application axe sufficient to readUy implement the method for multi-dimensional data 
compression as well. 

The included implementations have not been optimized. (Fbr instance, straight-forward 
optimization has improved run-time by as much es 70% over the included awft code; 
on the other hand, optimization destroys intelligibility of the codmg m relation to 
the underlying algorithm, and this explains the inclusion of unoptimized source code.) 
Another substantial improvement in run-time over the included source codes can be 
achieved by factoring so as to replace various manipulations of complex numbers with 
manipulations of integers as discussed before. 

Moreover, for clarity, the included source codes do not take advantage of renormal- 
ization. but for completeness, to each of awft and awift Is appended source code to 
replace the subroutines tgener and bgener by corresponding subroutines which do em- 
ploy renormalization. 
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